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ABSTRACT 

Firefly Algorithm (FA) is a newly proposed 
computation technique with inherent parallelism, 
capable for local as well as global search, meta¬ 
heuristic and robust in computing process. In this 
paper, Firefly Algorithm for Dynamic System 
(FADS) is a proposed system to find instantaneous 
behavior of the dynamic system within a single 
framework based on the idealized behavior of the 
flashing characteristics of fireflies. Dynamic system 
where flows of mass and / or energy is cause of 
dynamicity is generally represented as a set of 
differential equations and Fourth Order Runge-Kutta 
(RK4) method is one of used tool for numerical 
measurement of instantaneous behaviours of dynamic 
system. In FADS, experimental results are 
demonstrating the existence of more accurate and 
effective RK4 technique for the study of dynamic 
system. 

Keywords: FADS, Firefly Flashing, First Order 
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Introduction 

S. X. Yang of Cambridge University in the year 2008 
has proposed new nature inspired, swarm intelligence 
based optimization method on the idealized behavior 
of the flashing characteristics of the common fireflies 
like Photinus [lJ.This algorithm assumes that each 
firefly is a solution of the optimization problem and 
quality of the solution is proportional to the intensity 
of glowing, and consequently brighter firefly attracts 
its partners which may explore the search space for 


the better solution. Yang has shown higher efficiency 
and effectiveness for different multi-dimensional, 
nonlinear and multimodal optimization problems [2], 
In the year 2010 he has tested & designed several 
benchmark functions where this algorithm is 
outperforming other existing meta-heuristic 
algorithms [3], Fukasik et ah, 2009 has shown how 
FA performs better for continuous constraints 
optimization problems [4], After the introduction this 
efficient computation technique it is found in 
literature that different applications^] have developed 
to solve various common NP-Hard problems like data 
clustering[8], bi-level optimization[7], non-linear 
optimization problem [9], load forecasting problem 
[10] and multi-objective job-shop scheduling problem 
[12] .Different improvements based chaos [11] and 
there is hybridization with support-vector machine 
(SVM) [10] have been also emerging as more 
efficient optimization technique. Yang has shown in 
Fevy Flight distribution is outperforming than 
Gaussian distribution in randomized characteristics of 
fireflies [13].Dynamic systems are generally 
represented as a system of differential equation as 
parameters are varying with time and interrelated 
properties [14-20] and to measure different 
instantaneous property values the numerical solution 
of such differential equations are more practical and 
different numerical analysis are performed for the 
same [16,17]. In this present work FA has been 
studied and used to know the instantaneous behaviors 
of such dynamic systems modeled as system first 
order differential equations [18]. We are looking for 
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numerical solutions, so we are focusing on numerical 
analysis based solutions and successfully proposed 
more accurate formula for the RK4which has been 
derived using Taylor Series[19] and embedded 
successfully in our proposed FADS. 

Firefly Algorithm 

Firefly Colony 

Firefly Algorithm is most general optimization 
algorithm which can handle multi-dimensional, linear, 
non-linear, multimodal, non-convex, non- 
differentiable and continuous constraints optimization 
problem. Yang 2010 has shown that the Particle 
Swarm Optimization (PSO) is a special case it [3, 8 ], 
FA has designed after inspiring on the unique short 
and rhythmic flashing pattern generated to attract 
mating partners and / or potential prey by the most of 
the fireflies’ species. The radiation of this flashing 

light obeys the mathematical relation: I oc 1 ■/ 2 where 

I is the intensity of flashing light at distancer, also 
due to absorption this intensity again decreases and 
combined effect limits the fireflies to communicate 
within few hundred meter of distance. Yang 
formulated this flashing behavior as an objective 
function and developed a new bio-inspired, swarm 
intelligence based meta-heuristic optimization 
algorithm. 

Modeling the Attraction of Firefly 

This flashing behavior of fireflies can be idealized and 
following three idealized rules for the firefly colony 
system were proposed by Yang 2010: 

Rule I: All fireflies are unisex, so that one firefly is 
attracted to other fireflies regardless of their sex 

Rule II: Attractiveness is proportional to their 
brightness, thus for any two flashing fireflies, the less 
bright one will move towards the brighter one. The 
attractiveness is proportional to the brightness and 
they both decrease as their distance increases. If no 
one is brighter than a particular firefly, it moves 
randomly. 

Rule III: The brightness or light intensity of a firefly 
is affected or determined by the landscape of the 
objective function to be optimized. 


Using these three ideal rules and natural absorption of 
the media it is possible to formulate the attractiveness 
of the firefly as follows: 

P O') = 77771 ~ Po e ~ Yr2 (Gaussian form) 

( 1 ) 

where y is the absorption co-efficient (or attenuation 
factor) of the media and /? 0 is the initial attractiveness. 
From the inverse square law of light intensity it is 
possible to write for practical media: 

7(r) = ~~^2 ~ h e ~ Yr2 (Gaussian form) 

(T 

Here one is added to the denominator to avoid the 
divisibility by zero and I s is the intensity at light 
source. 

Modeling the Movement of Firefly System 

The Cartesian Distance between two fireflies Fj 
and Fj, are presently lying in n-dimensional search 
space at and Xj respectively is given in following 
equation (3) and used in the new position calculation 
formula: 

r ij = ||Xi - Xj || = V2L- */,fc) 2 (3) 

Fireflies are moving due to attractiveness generated 
for rhythmic pattern of flashing can be written as: 

X i+1 =X t + p 0 e- Yr ti{Xj - X t ) + a (rand - 7 ) (4) 

Here firefly Fj is the brighter or more attractive than 
firefly Fj, so firefly F t is moving towards Fj and the 
position of Fj i.e. X t is updating using equation (4) 
and it is also updating the intensity of Fj using 
equation ( 2 ) and hence the objective functions of the 
optimization problem. Attractiveness is updated using 
formula (1). For maximizing problems, intensity is 
directly proportional to the objective function value 
(/(A) oc f(X )) and for minimization it will be 
inversely proportional ( I(X ) oc l/f(X)). In this 
position changing equation second part is due to 
attractiveness while the third part is for 
randomization, like mutation in GA [23] to avoid 
stuck-at-local optima problem, a is the randomization 
parameter and rand is a random number generator 
uniformly distributed in [ 0 , 1 ]. 
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Algorithm for Simulation of Firefly Colony attractiveness with certain randomization due 

environmental interactions here is the pseudo code of 

Based on above mentioned modeling of attractiveness the firefly system to solve the optimization problem, 
and then the movement of firefly due to this 

Algorithm Firefly Algorithm (FA) 

( *maxf(X ) is the objective function,X = (x 1 ,x 2 ,..., x n ), position or solution vector in the constraint n — 
dimensional search space*) 

Step 1: [ Initialize ] 

Create N number of Fireflies F = (F 1; F 2 ,..., F ) 

N 

Place them at random positions S = (X t ,X 2 ,, X N )in the feasible solution region. 

Evaluate intensities / = (f, 1 2 ,..., I N ) for all N Fireflies 
Initialize light absorption co-efficient (y) 

Set maximum number of generation ( Gen Max ) 

Set current generation number (t = 1) 

Step 2: [Observe the swarming effect of the fireflies ’ colony for optimization ] 

while(t < Gen Max ) do 

begin 

fori = ItoAido 

begin 

for; = ltoido 

begin 

if(/ ; - < /;)then 
begin 

2.1 Move Fj towards F t in n — dimensional space using (3) 

2.2 Update light intensity using and Evaluate new solutions (2) 
end 

Update attractiveness using (1) 
end 
end 

Rank the fireflies using Sort operation 
Find the current best firefly 
Move all fireflies randomly in the search space 
Update generation counter t = t + 1 
end 

Step 3: [ Finished] 

Return best Solution Vector and Optimized Value 

Here the attractiveness depends on both flashing intensity and objective function value and distance among the 
fireflies which monotonically decreases with distance. But, for FA visibility is adjustable and more versatile in 
attractiveness variations, and this leads to higher mobility so the search space is explored more rapidly and 
execution time decreases. 


Dynamic System 

Dynamical system is a system in which any 
measuring parameter is varying with time. The world 
is full-of dynamic systems; all populations, human 
and otherwise, are dynamic systems, epidemics are 
dynamic systems, economies at all scales are dynamic 


systems, weather is complex dynamic system. 
Anywhere there are flows of mass and / or energy 
observed is a dynamic system. For the mathematical 
study of such varying system, instantaneous rate of 
changes parameters are modeled as system of time 
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dependent Differential Equations. To measure any 
such variable property or characteristic we need to 
solve all such differential equations after imposing 
boundary conditions. In all cases solving such 
equations are not so simple and / or in some cases it is 
impossible and generally numerical approximate 
solutions are preferred over complex symbolic 
mathematical expression [16-20],Euler, Euler- 


Cauchy, Taylor Series, Runge-Kutta Second Order & 
Fourth Order etc. methods most common to use for 
numerical solution. Runge-Kutta Fourth Order 
Method is most popular method because it is 
advanced form of general Taylor Series Method [18- 
20] and mostly used for the dynamic systems which 
can be modeled as first order differential equations. 


Example of Dynamic System 


Fluid system in nature is very common dynamic system, in the year 1872, Joseph Boussinesq[21,22] has 
proposed very simplified model for water waves known as solitary wave or soliton still now it considered as 
reference and valid for weakly non-linear and fairly long waves. Following system of first order differential 
equations [equations no. (5) & (6)] known as the Boussinesq equations arerepresenting this model. 


(ID 

dt 

dB 

dt 


= — e ~ 3 x B 

Pb 

= -D x B 

2 


D = {D 1 ,D 2 ,D 3 ) T 


(5) 

( 6 ) 

- is horizontal displacement andf? = (B 1 ,B 2 ,B 3 ) T is vertical displacement, 



0 

-1 0] 

dimensional constant and e 3 = 

1 

0 

0 


-0 

0 

0- 


—is a non- 

Pb 


d£i _ 

dt p b 2 


dD 2 _ g_g 
dt p b 1 



(?) 

( 8 ) 

(9) 


^ = \{D 2 B 3 -D 3 B 2 ) (10) 

^= 1 -(D 3 B 1 -D 1 B 3 ) (11) 

^ = i(£» 1 B 2 -D 2 5 1 ) (12) 

In above differential equations [equation no. (7) - (12)] we are showing component of different velocities.For 
time dependent instantaneous behavior these equations must be solved and use of RK4 has been shown by 
Desale et al in the year 2013 [21]. 


Runge-Kutta Methodfor Dynamic Systems 


To study the dynamic behavior of any system the gradient or variation relations of the system are expressed as a 
set of differential equation of equal or different orders. Most of the simple form is the first order form called as 
Ordinary Differential Equations (ODE). The solutions of these ODEs describe the functional behavior of the 
system and for this purpose different calculus oriented methods like integrating factor etc are used. The RK4 
method is very efficient, as it avoids all kind of differentiations required in Taylor Series method. In the 
derivation of the RK4 method we get eleven independent algebraic relations for thirteen different coefficients 
and these are forming a system of equations with rank less than number unknowns. To find the value of 
unknowns Runge-Kutta had assigned empirically fixed values for two such variables and then solve the system 
of equation using numeric algorithm [19], finally we got the most popular recursive equation for evaluation.Let 
for a dynamic system there are N numbers of properties to measure and these are represented by the time 
dependent vectorT(t) = {y l (.t)}i =1 and F 0 = Y (t 0 ) = {y f (t 0 )}[f 1 be initially measures or know parameters. 
Then the instantaneous dynamic behavior represented by 1 st order differential equation form as follows: 
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= F{t, Y(t)) = V(t))}fL 1 with initial condition (t 0 , T 0 ) (13) 

And the solution will be of the formF = G(t) = Different complex forms of F(t, Y (t)) are possible 

and to get G(t) different techniques are used and for numeric solution different approximate numeric 
algorithms are used. Most popular numeric method is RK4. This method is derived from Taylor Series method 
after approximating it up to 4th order term in the series and Runge-Kutta had assumed the recursive solution of 
the following form [19, 20 ]: 


y l k+ 1 = y k + w ± k[ + w 2 k l 2 + w 3 k 3 + w 4 k l 4 , where (14) 

K = hfi(t k , y k (t), y k (t), ...,y k (t) ) (15) 


k l 2 = hfi(t k + a-tKyl + b^lyl + b x k\, + b-^kf) 


(16) 


H = hfi(t k + a 2 h,y k + b 2 k[ + b 3 k[,yl + b 2 k\ + b 3 k\, ...,y% + b x k\ + b 3 k l 2 ) (17) 


k\ = hf{t k + a 3 h,y k + b 4 k[ + b s k l 2 + b 6 k l 3 ,y% + b 4 k[ + b s k\ + b 6 kl ...,y% + b 4 k[ + b s k[ + b 6 k l 3 ),Vi = 
1. N (18) 


After coefficient comparison with the Taylor’s series method following eleven equations for thirteen 
unknowns{ai)f =1 , { 6 ;}f =1 and{Wj}f =1 of above equations (14) - (18) have generated, for which any conventional 


method is not suitable to solve. 

b 4 = a ± , b 2 + b 3 = a 2 , b 4 + b 5 + b 6 = a 3 (19)-(21) 

1 

w 1 + w 2 + w 3 + w 4 = l,w 2 a! + w 3 a 2 + w 4 a 3 = - (22)-(23) 

w 2 af +w 3 al + w 4 a 3 =-,w 2 cl\ + w 3 a| + w 4 a 3 = - (23)-(25) 

3 4 

w 3 a 4 b 3 + w 4 (a 1 h 5 + a 2 b 6 ) = \,w 3 a 1 a 2 b 3 + w 4 a 3 (a 1 h 5 + a 2 b 6 ) = \ (26)-(27) 

D O 

w 3 a\b 3 + w 4 (alb 5 + ajb 6 ) = w 4 a ± b 3 b 6 = (28)-(29) 


These eleven equations [equations no (19)-(29)]are forming a system of equations for which we have thirteen 
unknowns and eleven relations, so the rank is less than number of unknowns [18]. Runga-Kutta equalized the 

rank after assuming^ = -,b 1 = 0 , then solved these equations and finally derived the popular recursive 
relation from equation (14) as: 


y l k+i =y l k + + 2k 2 + 2k 3 + k 4 ) (30) 

11 11 
with the solutionsa 2 = a 3 = 1 , b 2 = b 3 = b 4 = b 5 = 0, b 6 = 1, w 1 = w 4 = w 2 = w 3 = - 

2 2 6 3 


Firefly Algorithm for Dynamic System (FADS) 

For the design of FADS for general dynamic system we use these thirteen parameters as 13-dimensional 
position vector: X = (a 1 ,a 2l a 3 ,b 1 ,b 2l ... ,b 6 ,w 1 ,w 2 ,w 3 ,w 4 ) E R 13 for the different fireflies of FA. Above 
equations (19) to (29) are nothing but equality constraints to satisfyfor best possible position vector (X) to get 
the more accurate recursive formula represented in equation (14). For fitness or objective function calculation 
thesquare difference of left hand sidealgebraic expressions of these thirteen parameters with right hand side 
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constant for these eleven equations for certain position vector or firefly is considered as square error 8 t and sum 
of these square errors as objective function and we have to minimize it towards zero value: 

MinimizeI(X ) = (31) 

For Algorithm Firefly Algorithm (FA) mentioned above we have incorporated position vector of the fireflies 
along with square error measurement in equation (31) as intensity of these fireflies. After getting more accurate 
13-parameters values recursive formula (14) with other formulas for the calculations of k/s mentioned 
equations (15)-(18) are used for numerical values of characteristics parameters on time for a dynamic system. 


Experimental Results 

Proposed FADS has been implemented using C++ programming language and tested with FA related 
parameters setting: Number of Firefly (N) = 20, Number of Generation ( Gen Max ) = 1000, /? 0 = 0.2 , and 
Y = 1.0. This has been executed for 30 independent runs with said parameter values and we get different 
values for these thirteen variables which were found by Runge-Kutta and these are shown in Table 1. From this 
we can conclude that here is a new form of the recurrence relation which is outperforming the existing system. 

Table 1. 


Parameters 

RK Solution 

Average FADS Solution 

a 4 

0 

0.197491 

a 2 

1 

2 

0.665014 

a 3 

l 

0.813665 

h 

1 

2 

0.201191 

b 2 

1 

2 

-0.737769 

^3 

1 

2 

1.409750 

b * 

0 

1.005670 

b* 

0 

-0.864207 

h 

l 

0.667915 

W 1 

1 

6 

0.109450 

vv 2 

1 

3 

0.277182 

w 3 

1 

3 

0.362196 

w 4 

1 

6 

0.248217 


To study the performance of this dynamic system analyzer we have tested an expontial type system defined by 
the following first order differential equation: 

~dP~ = = ^(0) = ty(t) with initial condition (t 0 , T 0 ). 



Fig. 1. 3D Comparative Study 
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The effectiveness of the system is presented by a comparative study of existing solution and standard C++ 
language library with the proposed FADS and this has been shown graphically in above Fig. 1. 



The convergence of FA to achieve the best possible parameters values is also shown in Fig. 2. 


Conclusion 

This work has proposed how the robust and efficient 
FA optimization technique can be used to get numeric 
solutions of complex dynamic system represented by 
a set of differential equations form. This also suggests 
more accurate alternative of RK4 and also suggesting 
simpler way to develop RK-N where N > 4 for more 
accurate behavioral study of complex dynamic 
system. 
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